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Abstract 

The authors provide the formulation basis for nonlinear transient 
analysis of finite element models of structures using energy minimiza- 
tion. With both geometric and material nonlinearities included, the 
development is restricted to simple one and two dimensional finite ele- 
ments which are regarded as being the basic elements for modeling full 
aircraft-like structures under crash conditions. The results presented 
establish the effectiveness of the technique as a viable tool for this 
purpose . 
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I. Introduction 


Nonlinear transient analysis of structures has been of increasing 
interest to engineers by virtue of their motivation for minimizing human 
and property damage resulting from the catastrophic failure of such 
structures under crash or seismic conditions. Complexities of the 
structural configuration and its equally complex transient response in 
the presence of material inelasticity make finite element modeling of 
such structures a very natural and plausible recourse. Portions of the 
structures may remain elastic and undergo infinitesimally small deforma- 
tions while other portions may experience finite deformations and mo- 
tions and respond inelastically under time-varying loads that may lead 
to a complete failure of the structure. If finite strains are to be 
permitted in the model, distinction must be made between undeformed and 
deformed configurations and the concepts of pseudo stresses and conjugate 
strain measures which have intricate physical interpretations must he 
introduced [1], Furthermore, strictly speaking most elastic-plastic 
theories which hypothesize an additive decomposition of the total strain 
into an elastic and a plastic component lose their validity in the large 
strain domain [2], Because of this, most developers of nonlinear analy- 
sis codes restrict themselves to a small strain formulation but permit 
finite displacements and rotations thereby ailowing buckling and collapse 
of the structure to occur. There are some indications that this may be 
adequate for most practical purposes. 

With this hypothesis as its basis, the present discussion focuses 
on the simulation of response of a structure modeled as an assemblage 
of membrane, frame (3-D beam), stringer elements and rigid links (see 
Figure 1). The mathematical model is a finite element displacement 



Figure i 






model which consists of discretizing the actual structure by an assem- 
blage of finite elements and approximating the response of each element 
by a finite number of deformation states expressed as linear functions 
of the generalized nodal displacements. 

Two distinct solution approaches exist: (i) the vector approach 

and (ii) the scalar approach. In the former , the mathematical model is 
derived on the basis of the principle of virtual work and reduces to a 
system of nonlinear second-order differential equations in time. In 
the latter approach , a scalar or potential function associated with the 
energy of the model is introduced, minimization of which yields the 
desired equilibrium configuration. In both approaches a temporal finite 
difference scheme is utilized to effectively eliminate time as a var- 
iable. As a result, in the vector approach the equations of motion are 
reduced to a system of nonlinear algebraic equations in the unknown 
nodal parameters of the finite element model^ In the scalar 

approach, which is of relevance to this paper, the problem Is reduced to 
a well known problem in mathematical programming namely the unconstrained 
minimization of a nonlinear function of several variables. 
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II, Minimization Technique for Nonlinear Analysis 


a. Formulation Basis 

The scalar approach has been successfully used for static and tran- 
sient nonlinear structural analysis [7]-[9]. In this case the problem 
of response prediction is equivalently posed aa the minimization of a 
potential function of the unknown nodal parameters of the finite ele- 
ment model. For all structural problems with geometric and material 
nonlinearities of the type considered herein such a potential function 
always exists. Although this technique has been hitherto used for 
mainly positive or negative definite systems, other systems which fail 
to be positive or negative definite can be handled by using the least 
squares method or the modified conjugate gradient method with precon- 
ditioning [10]. In some cases for such systems displacement incremen- 
tation rather than load incrementation in conjunction with conventional 
unconstrained minimization techniques can also be equally effective 
111 ]. 

The minimization scheme as applied to the solution of transient 
nonlinear structural analysis problems consists of minimizing a poten- 
tial function associated with the system for an assumed relationship 
between displacements and time. The displacement-time relation for each 
generalized nodal displacement of a finite element model may be assumed 
of the form [12] 

- 6(At) 2 q el + <i - e)(At) 2 q 01 + <At)q 01 + q 01 (1-a) 

q ei - Y(At)q el + (1 - Y)(At)q 0i + q Qi (1-b) 

where q fi ^ is the i-th generalized nodal displacement at the end of the 
time step and 0 and y are constants. These constants are determined in 
terms of the i-th generalized nodal displacement, velocity, 

and acceleration, <jg at the beginning of the time step and the generalized 
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nodal displacement, q^, at the end of the time step. It can be easily 
verified that the equation of equilibrium for an N degree of freedom 
system with lumped masses 


M l%i - F l + 


; i - 1,2. . .N 


correspond to the necessary conditions for the functional 


i?i (l 26Wt) 2 " ( B(L) 2 ,0i + S <*‘> ^ + <2S 


- F i (t+At)q gi } +U+C (3) 

to be stationary. In Equation (3), U is the strain energy and C is an 
arbitrary constant. Thus, knowing q Qi , q Qi and q Qi at time t for any 
given load at time (t+At), the functional S may be minimized with 
respect to the generalized nodal displacements, q gi (i-l,...N), in order 
to determine the corresponding stable equilibrium configuration. Thus, 
this scheme satisfies equilibrium at the end of the time step, thereby 
providing an implicit temporal integration scheme. The size of the time 
step is automatically controlled so that the error at half time based 
on interpolated configuration data is less than a prescribed change in 
total energy. In general, the strain energy U will be a nonlinear 
function (at the very least a quadratic) of the generalized nodal dis- 
placements q g ^. Details on the explicit evaluation of U as a function 
of will be touched upon later. 

Of all the available techniques for unconstrained minimization only 
tiie quasi-Newton or the variable metric methods have been more frequently 
used for structural analysis, because of their higher effectiveness 
[13]. Again, unless one accounts for the sparsity of the variable 
metric, one has to almost invariably resort to some form of a conjugate 
gradient technique for problems wherein N Is an extremely large number. 
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The extension of the minimization techniques to extremely large scale 
nonlinear structural analysis problems is a subject of separate research 
[14] in itself and is beyond the scope of this papi . 

Most algorithms for unconstrained minimization seek a direction of 
travel and the amount of travel in that direction. The manner in which 
theBe are sought depends upon the sophistication of the particular 
algorithm invoked. Most often the directions of travel are sought in a 
manner which guarantees not only a decrease in the value of the function 
to be minimised at each iteration but also a convergence to the minimum 
in a finite number ot iterations (usually N+L for an N dimensional 
space) in the case of quadratic functionals. It is important to note 
that ail functionals are very nearly quadratic in the nei ghborho- ><i of 
the minimum. The iterative scheme is begun with an initial guess which 
is usually the null vector in the absence of other better estimates. 

For the variable metric or the conjugate gradient methods the required 
gradient of S is evaluated either analytically or by a finite difference 
operation on S. The use of an analytic gradient results in a substan- 
tial saving in computational effort. This saving is the result of not 
only a cheaper gradient evaluation but most often a faster convergence 
of the solution because of higher accuracy of all computed quantities 
[13]. The i-th component of the gradient of S can be written as 


3S M •• p , 3U 
tt * M.q . - F.+ -5 

a, el 1*1 1 3 q ei 

The term in Eq. (4) requiring significant computational effori is 


(4) 


— as it embraces the geometric and materia] nonlinearities. Using 

half-station central differences this is given by 
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* • * 


3U 

3q 


V.r q u 2 ' 


<*.1 + 2 %!+! 


q eN* 


el 


Aq 


•i 


f U k (q el' q e2 
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q ei ‘ 2 Aq ei* q ei+l *“ q eN ) 


Aq 


( 5 ) 


ei 


where Aq g ^ is a small change in the i-th component and m is the number 
of members or elements which has the i-th degree of freedom in common. 
In evaluating the gradient vector analytically [17], each of its com- 
ponent involves the evaluation of only a single function similar to the 
function for member energy evaluation. Thus, 

3U f 9W 


3q 


ei 


ft/ aiT d V ft / 


,dWv / 3e . . 

( “ ) k ( 3T” ) k dv k 
k aq k k 


(6-a) 


where 


W ■ strain energy density 

2 


73 (eL + E. 2 .. “ E.„e„.. + - yL ) 1/2 for a 


xx 


£ ■ effective strain 


yy xx yy 4 xy 

two dimensional stress state 


( 6 -b) 


« £ for a uniaxial stress state 
xx 


and for one step incremental loading or unloading 


<^> 

d£ 


effective stress 


- la 2 + o 2 - a a + 3T 2 ) 1/2 

xx yy xx yy xy 

for a two dimensional stress state 

* o for a uniaxial stress state 
xx 


< 6 -c) 


Equations ( 6 -a) through ( 6 -c) imply the use of Hencky's total strain 
theory along with its assumption that fn the strain hardening r mge che 
inelastic component of the total strain is predominant [15], This in a 
way is consistent with the assumption that the total strain can be de- 
composed into an elastic and a plastic part especially in cases where 
the strains are large [2], According to reference (2] it is only when 
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plastic strains are predominant that such a decomposition Is justified. 

The problem at hand could have equally well been formulated using 
the Incremental flow theories of plasticity in the strain hardening 
range. The potential function instead of being a function of the total 
quantities need then be expressed in terms of incremental quant ies 
and the minimization technique can still be used [16]. As a matter 
of fact, it may be conjectured that the performance of the solution 
algorithm will perhaps be significantly improved using such a formula- 
tion even though the material model may then be slightly more complex. 
In any event, it is immediately obvious that a significant reduction in 
computational time will be realized if analytic gradients are used in 
preference to central difference gradients. 

The complexity of the strain energy evaluation for any clement is 
determined by its deformation model. This is discussed next, 
b. Deformation model : 

The deformation model of the entire structure is synthesized from 
deformation states of each element of the structure. These states are 
expressed in terms of generalized displacements of the nodes of the 
structure at which the elements interface. 

The displacement field within each element is chosen as a continu- 
ously differentiable function of the local spatial coordinates and the 
generalized nodal displacements. The field maintains intcrelement 
continuity of its essential derivatives thereby providing a Galorkin 
model of the system. The local generalized nodal displacements of each 
element are then related to the global displacements of the assemblage. 
These relations, which can be interpreted as transformations of the 
local coordinate system to the global coordinate system, may be linear 



or nonlinear depending upon whether the motions and deformations of the 
elements are infinitesimal or finite. For large rigid body rotations, 
these transformations are accomplished using Ruler angles which are 
linearly independent by virtue of the fact that the rotations are per- 
formed in a prescribed order. 

There are three kinematic descriptions most commonly used for char- 
acterizing large displacements of finite element models of structures. 
These are: (i) the total Lagrangian formulation wherein the initial 

undeformed configuration is the reference configuration, (ii) the up- 
dated Lagrangian formulation which uses a total Lagrangian formulation 
within each load or time step but updates the reference conf iguraticn at 
the end of each step and (ill) the co-rotational or rigid converted 
coordinate formulation which utilizes a coordinate system rigidly at- 
tached to an element and moving with the element. For development of a 
large rotation formulation vital for crashworthiness studies the use of 
the total Lagrangian formulation is unsuitable since most structural 
theories permit only moderately small rotations [18], The co-rotational 
formulation decomposes the total displacements into a rigid body motion 
component and a strain producing component. Thus, with the restriction 
of small relative rotations within the element, this formulation leads 
to a simplification of the strain-displacement relationship on the 
element level while still permitting arbitrarily large rotations of the 
element. The present deformation model uses the co-rotational or rigid- 
eonvected formulation for its kinematic description. 

Through appropriate kinematic constraints modeling of massless 
degrees of freedom or of deformat ion- free rigid links or even the simu- 
lation ot contact with an impenetrable, rough plane are easily achieved. 
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Rigid link* can be used to simulate either joint eccentrices or rigid 
parts of a structure. In the interest of a truly unconstrained minimi- 
zation Lagrange multipliers or penalty functions are avoided. Rather 
kinematic constraints are formulated as prescribed displacements un-i. r 
reactive forces provided by the gradients of the strain energy with 
respect to the corresponding degrees of freedom, 
c. Material model : 

Although closed form analytic expressions for U can be developed 
when the material is elastic the same is not true when elements yield. 
Then the response depends upon the current values of stress components 
and the past history. Von Mices' yield criterion together with Hem-key 's 
total strain theory provides a simple means of calculating strain 
energy density distributions throughout an element that has yielded. 
Because total stresses and total strains are no longer linearly related 
recourse must be made to numerical integration (Gaussian or Lobatto) of 
the strain energy density over the volume of the element. Thus, a frame 
element which was strictly a uniaxial member in the elastic range, 
typified by its cross-sectional area and moments of inertia, requires a 
full three dimensional characterization in the inol > tic range. In 
other words, frame elements for inelastic analysis require a classifi- 
cation based on the different cross-sections. This development is 
restricted to frame elements with thin-walled sections of the closed and 
open (Box, Tube, Elip and E) variety - a characteristic of general 
aviation aircraft frames. However, in the elastic range the development 
does permit frame elements with arbitrary cross-sections characterized 
by their gross section properties. In the interest of simplicity, 
classical shear flow theory for thin-walled sections is used and certain 
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simplifying assumptions regarding torsion, warping and shear deforma- 
tions in the inelastic range are introduced. This is characteristic of 
most nonlinear analyzers mainly because the development of a truly three- 
dimensional frame element for nonlinear inelastic response is a formid- 
able task perhaps even more challenging than that of the development of 
a plate bending or a shell element for the same purpose. In fact, in 
the inelastic range it may be easier to model a thin-walled beam of 
arbitrary cross-section by an assemblage of a large number of plate and 
shell elements thereby permitting a faithful representation of very 
complex effects like restrained warping, torsion, cross-sectional dis- 
tortions, etc. 

Thus, it is clear that when an element yields, the complexity of 
the strain energy evaluation Increases several times in relation to its 
purely elastic behavior. A number of quadrature points have to be 
assigned over the volume of the element and using the material model 
stresses and strain energy densities have to be evaluated at each of 
these points for known values of strains (Figure 2). The average strain 
energy density which is simply the weighted sum of these strain energy 
densities then enables the calculation of the total strain energy. It 
must be noted, however, that the stress-strain history at each of these 
quadrature points, which corresponds to a unique location on an ideal- 
ized effective stress-effective strain curve for the material of the 
element (Figure 3), must be made available at all times. This places 
highly increased demands on computer storage as inelastic deformations 
progress with time. 

The material is assumed to unload elastically. For modeling plas- 
ticicy under cyclic loading kinematic hardening with an idealized 
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Figure 2 
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Bauschinger effect is assumed. Specialized elements like the gap ele- 
ments and stays can be easily modeled simply through an appropriate 
modification of the material model of the conventional elements. 

The details of the total strain energy calculations for the dif- 
ferent element types considered and the transformations relating element 
behavior to global variables is relegated to the appendix. 
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III. Results and Discussion 


The effectiveness of the minimization technique in solving non- 
linear problems is very much a function of not only the size of the load 
or time step but also the extent and type of the nonlinearity-geometric 
or material and even the type of the temporal discretization scheme used 
which is to say the assumed values of $ and y in Eq. (1). With this in 
mind, the process of selection of problems for validation was geared 
towards providing an evaluation of the techniques under different types 
of nonlinearities. Problems belonging to three distinct classes namely: 
(i) quasi-static, elastic with geometric nonlinearities, (ii) quasi- 
static, elastic-plastic with geometric nonlinearities and ( iii ) transient, 
elastic-plastic with geometric nonlinearities were selected. Independent 
solutions or experimental results for these problems were available for 
comparison purposes. 

Figure 4 shows the case of a rod-spring problem wherein the stiff- 
ness of the spring is just enough to prevent a snap-through and provide 
a single-valued load deflection response. Most researchers regard this 
problem as geometrically highly nonlinear. Using stringer elements with 
load steps as high as 1 lb. , the energy minimization solution is indis- 
tinguishable from the easily obtainable exact solution to this problem. 
Higher load steps could have been chosen but caution must be exercised 
with extremely large load steps since the performance (the number of 
minimizations required for convergence) of the minimization algorithm 
may be adversely affected. In other words, the computational effort 
within a load step may increase substantially enough to offset the sav- 
ings accured from fewer load steps. 
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VERTICAL DISPLACEMENT IN 


Figure 5 provides yet another example of such a nonlinearity except 
that in this case the load response curve is no longer single valued but 
is a composite of stable and unstable branches. Using straight-forward 
load incrementation, it is possible to locate only the stable equili- 
brium configurations as Indicated in Figure 5. Using displacement 
incrementation, however, the entire load response curve can be easily 
obtained. The response predicted by energy minimization agrees very 
closely with that predicted by the nonlinear analyzer developed by 
Stricklin and Haisler [19] for an identical model of the shallow arch 
using frame elements. 

While both of the previous problems involved only geometric non- 
linearities, Figure 6 presents the case wherein both material and geo- 
metric non Linearities interact. The experimental prediction of the 
post-buckling, elastic-plastic response of this beam-column with a thin- 
walled channel cross-section was the result of a test carried out by 
Anderson et al [20], To prevent failure by direct compression, the 
column was tested at an inclination of 5° from the vertical with both 
ends of the column being clamped. A column under these conditions is 
highly imperfection sensitive and hence Anderson et al assume an addi- 
tional 1° offset, as shown in the figure, for the mathematical model 
hoping to simulate the inherent imperfections of the actual column 
tested. Because the response involves a highly unstable branch, dis- 
placement incrementation had to be used in place of load incrementation. 
The response predicted by energy minimization agrees extremely well with 
the experimental prediction and even more so by comparison with that 
predicted using UMVCS-1[21]. 
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VERTICAL DEFLECTION 



Figure 7 illustrates the case of the transient response in the 
presence of both geometric and material nonlinearities. The impulse is 
large enough to cause the entire beam to respond inelastically while 
experiencing moderately large relative rotations. The experimental re- 
sponse for this beam was obtained by Krieg et al [22]. It is immediately 
obvious that the quality of the response prediction is very much a func- 
tion of the values of 6 and y. This is not to say that optimum values 
of 0 and y «*xist which guarantee optimum fidelity of the response pre- 
diction. Ail a matter of fact, optimum values of 6 and y appear to be 
very much problem dependent. Again, the response may be significantly 
affected by the use of a consistent mass matrix and with rotatory iner- 
tia and shear deformation effects included. 

With the possible exception of the problem of Figure 8, all the 
previous problems involved only a relatively few degrees of freedom. 
Furthermore, the state of stress in all of these proble is was essen- 
tially uniaxial for all practical purposes. Using constant strain 
membrane elements the maximum strain in the vicinity of a notch in the 
direction of loading is determined and compared with the experimental 
results. The agreement between the two predictions is good but could 
perhaps be improved upon by the use of nonlinear strain displacement 
relationships in the co-rotational coordinate system. 

This demonstration of the effectiveness of the minimization tech- 
nique as a tool for nonlinear analysis has been, no doubt, restricted to 
classical, idealized problems with a relatively few degrees of freedom. 
Extensions to large scale problems like the section of an aircraft, as 
in Figure 1, or even a full aircraft may involve several thousands of 
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degree* of freedom. The effectiveness of the present technique for 
response prediction of such structures remains to be demonstrated. 

Using preconditioned conjugate gradient technique or variable metric 
methods which exploit sparsity, it is believed that this is no longer an 
insurmoun table task. Its performance vis-a-vis the pseudo force and 
hybrid techniques for such large scale problems will be the subject of a 
follow-on paper. Incldently, however, for small scale problems of the 
type considered herein the energy minimization technique has been shown 
to be at least comparable to, if not better than, the pseudo force 
technique [14]. 
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Appendix 

Evaluation of the Total Strain Energy 
The scalar approach for the solution of problems of structural 
analysis requires that the strain energy of the system be expressed, 
explicitly or implicitly, as a function of the global generalized nodal 
displacements of the finite element model. 

From a known vector of the generalized nodal variables in the 
global co-ordinate system, consistent with the prescribed boundary 
conditions, a vector of local generalized variables in the co-rotational 
co-ordinate system of each element is established through transforma- 
tions which are functions of its geometry and its rigid body rotations. 
The assumption of deformation patterns of the element as functions of 
these local generalized nodal variables (interpolating polynomials) 
yields element strains. Recourse to the element material model then 
yields the corresponding stresses and strain energy densities at various 
predetermined points (quadrature points) over the extent of the element. 
Barring purely elastic response, a simple weighted summation of thes^ 
quantities over the element volume yields stress resultants and strain 
energies respectively. For purely elastic response these are provided 
by well-known closed form expressions. For the elastic-plastic response 
the strain energy density may be decomposed into an elastic part and an 
incremental-dissipative part thereby providing an estimate of the total 
energy of the system that has been dissipated through inelastic deforma- 
tions. Thus, as shown in Figure 3 for a system with M elements 
M . M . M r 

U - £ U - £ U 1 + AUt - Jj ( J W £ dv + J AW d dv) (A-l) 

i-1 i**l 6 1=1 v i v i U 

where the dissipative energy AU^ is the incremental dissipative energy 
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computed from the previous stress state typied by the point a Q on the 
effective stress-effective strain curve of Figure 3. 

In the following sections expressions for the strain-displacement 
relations are developed for the different elements. 

(i) Stringer Element 

A structural component of uniform cross section which is initially 
straight and which is capable of resisting only axial loads is known as 
a stringer element. 

From Figure (A-l) it can be seen that for assumed nodal displace- 
ments (11^, V p , Wp) and (U^, V^, W^) of nodes p and q, the change in 
length, DL, of the element is given by 

DL - [(X q + U q - X p - U p ) 2 + (Y q + V q - Y p - V p ) 2 

+ < Z , + M , - 2 p ' V 2ll/2 - UX q - V 2 + % - Y P >2 

+ (z q - Z p ) 2 J 1/2 (A— 2 ) 

which can be simplified to 

DL - L[1 + + g . W , ,1/2 . J (A . 3) 

L L 

A being the difference operator for q and p end values. Assumption of 
the usual linear interpolation function in the co-rotational co-ordinate 
system then yields 

£ - (~) (~) (A-4) 

(ii) Frame Element (3D Beam) 

A frame element (Figure A-2) is a structural component which Is 
initially straight and which undergoes axial, bending and torsional 
deformations resulting from finite displacements and rotations of its 
ends. 


28 




From Figure (A-3), the displacements of the end q relative to the 


end p can be seen to be 


6u - (R - R ) - L + (U - 1J ) 

- -q -p -q -p 

or in terms of the three components as 


(A-5) 


1 6u 

| 

(X - X 

i q p 

I j L j 


U - U 
1 q P 

\6v 


Y - Y 

| q p 

" J°| 

+ tn p 

V - V 

| q p 

'6w 

! 

Z - 2 

A n ' 

( o , 




(A-6) 

where again U^, and (i*p or q) denote the global displacements of 

the nodes. The matrix [T] can be shown to be [9] 

P 

UHT,(9«. 0 . e )] 


[T]. 


l' T x’ 


with 


'VWV 1 


c c 

y z 


1 *P yp’ zp' 


c s 

y z 


(A-7) 


-C S +S 8 C c C +S S S 
xz xyz xz xyz 


-s *i 

y 

S C 
X y 


Ls s +c s C 
x z xyz 


-S C +C S S C c 

x z xyz x y 


(A-8) 


c - cos at. and s, * sin a. for i=x,y and z. Angles $ , <p and <p are 
i li l xyz 

the initial orientation angles described in Figure (A-2) and angles 

9 , 9 and 9 are the rigid body rotations of the end p. In deriving 

xp* yp zp 

Eq. (A-7) Euler angle transformations are implied with the order of the 

rotations being a , a and a . 

° z y x 

Similarly, with the restriction of small relative rotations within 

the element, the rotations \p , and ip of the end q relative to the 

x y z 

end p are 


(A-9) 


With the relative generalized displacements {<5u, 6v, 6w) and 

{ip , \p , tp } known the usual deformation patterns of the reference axis 
xyz 

of the beam element in the co-rotational co-ordinate system are assumed 


(ip 1 


9 - 0 \ 

x i 


1 xq xp 1 

\p 

1 3 [T] 

|9 - 0 ) 

y | 

P 

| yq ypi 

K' 

1 

>0 - 0 ) 
in 9n 


.11 



Figure A-3 
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Co be 


v(U - -J (3{ 2 - 25 3 )(Sv - z ji > + U 3 - $ 2 )«, I 

* I (A-10) 

w<£) - £ <3t2 ' 2 ? 3 >< Sh + >’A ) ‘ (?3 ‘ ( 

8 * »« ' 

where £ ■ x/L and y and z are the co-ordinates of the shear center of 

s s 

the cross-section of the beam. The strain of the reference axis can 
then be shown co be 

e - C % - nfrr (l-20(6v-*J> > + 2(3£-l)i|M 

\ 8 X Z (A-ll) 

- Clf (i-^)(6w+y^ x ) - 2<35-W y ] 

with rj ® y/L and C “ z/L. In the above equations it is implicitly as- 
sumed that the lateral displacements and twists are referenced to a 
longitudinal axis through the shear center while the axial displacements 
and rotations are referenced to the centroidal axis. As shown in ref- 
erence [23] this assumption necessitates the introduction of an addi- 
tional degree of freedom in the axial direction in the interest of equi- 
librium satisfaction in the inelastic range. 

(iii) Membrane Element 

The membrane element of Figure (A— 4 ) is a plane triangular thin 
plate element under constant strain. The element can undergo large 
rigid body motions but its deformation is restricted to only in-plane 
stretching resulting from finite displacements of its vertices. 

The orientation of the element is uniquely determined by the global 
co-ordinates of its three vertices, p°, q° and r°. The relative dis- 
placements <5u , 5u r and 6v^ defined in Figure (A-4) can be seen to be 
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C A— 12 ) 


<5u - R - R* 

q 

6u r ■ Q cos B - Q* cos a 
6u r - Q sin B - Q* sin a 


with 


R* - (X 2 + Y 2 + Z 2 ) lf2 

' rp rp rp 

* - « X rp + V* + % + V rp> 2 + < 7 ‘rp + V'’ 1 ' 2 
Q° - (X 2 + Y 2 + Z 2 ) l/1 

qp qp qp 


q - [(x + u ) 2 + (Y + v ) 2 + (z + w ) 2 ] 1/2 

qp qp qp qp qp qp 


(A- 13 ) 


cos a - (X X + Y Y + Z Z rn )/(Q"R°) 
qp rp qp rp qp rp 


cos B - t(X rn + U )<X rn + U ) + (Y + V ) (Y + V ) + 
r .P qp rp rp qp qp rp rp 


(z + W )(Z + W )]/(QR) 

qp qp rp rp 


and typically ot^ ■ “ a j* Next, as in the case of the two previous 

elements, deformation patterns u(x,y) and v(x,y) in the eo-rotational 
co-ordinate system when expressed in terms of the local nodal displace- 
ments yield 

6u 6u 6u 

u(x,y) " u p + ( )x + I cotajy 


6v <5v <5v 

r(x,y) - v p + ( -^ )x + [ cota]y 


(A-14) 


with the strains e , E and y determined on the basis of the small 
xx yy xy 

deformation theory as 

~ 6u 

£ » la . ( s \ 

e xx 3x { R° ' 


3v ^ V r 


yy 3y Q°sina 

* - 6u 6u 

y . 1h + ^ a cota) 

Y xy 3y 3x { Q°sina R* ota '’ 


(A-15) 
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Thus, with the assumption that the total deformation theory of plas- 
ticity is applicable, the effective strain and effective stress defined 
by Eqa. (6— b) and (6-c) yield estimates of the stresses and strain 
energy densities from the material model. It is obvious that the inte- 
grations over the volume of the element are rendered trivial by virtue 
of the assumption that strains and hence the stresses and strain energy 
densities within the element are constant. 

(iv) Rigid Link 

Rigid link is an element which merely translates and rotates 
without any appreciable deformations. The element is identified by two 
nodes located with reference to the global co-ordinate system. One of 
these two nodes is referred to as the master or primary node, p, with a 
maximum of six independent degrees of freedom. The motions of the slave 
or secondary node or nodes, q, are determined purely from kinematics by 
setting the left hand side of equation (A-5) to zero. Knowing the 
dependent displacements of the secondary nodes, 

{u} - {u} + [T]*{L} + {R> - {R} (A-16) 

q p p p q 

the contribution, to the total potential energy, of loads applied 
directly at the secondary nodes can thus be determined. 
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Figure Legend 
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